R Notebook

This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.

Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Ctrl+Shift+Enter.

library(spatstat)
library(sp)
library(rgeos)
library(maptools)
library(GISTools)
install.packages("tidyverse")
library(tmap)
library(sf)
library(geojsonio)
library(tmaptools)
library(tidyverse)
library(rgdal)
library(raster)
library(SciViews)
library(base)

####Basemap
library(tmap)
BoroughMap <- geojson_read("https://raw.githubusercontent.com/ft-interactive/geo-data/master/uk/london/london-wards-2014.geojson", what = "sp")
BNG = "+init=epsg:27700"
BoroughMapBNG <- spTransform(BoroughMap,BNG)
WGS = "+init=epsg:4326"
tmap_mode("view")
tm_shape(BoroughMapBNG) +tm_polygons(col = NA, alpha = 0.5)

#aggregate city of london 
library(sf)
library(dplyr)
library(rgdal)
BoroughMapSF <- st_as_sf(BoroughMapBNG)
library(rgeos)
CityofLondon<- gUnionCascaded(BoroughMapBNG[c(630:654), ])
plot(CityofLondon)
CityofLondonSF <- st_as_sf(CityofLondon)
CityofLondonSF$gss_code_ward<-"E05000001"  ##add same columns in city of london, for aggregating with other borough
CityofLondonSF$gss_code_borough <- "E09000001"
CityofLondonSF$borough <-"City Of London"
CityofLondonSF$ward <- "City Of London"
BoroughMapSFnew <- BoroughMapSF[-c(630:654), ]
LondonWardSF <- rbind(BoroughMapSFnew,CityofLondonSF) ##Attention New not new
LondonWard <- as(LondonWardSF,"Spatial")
rownames(LondonWardSF) <- 1:nrow(LondonWardSF)
plot(LondonWard)
tm_shape(LondonWard) +tm_polygons(col = NA, alpha = 0.5)

####population density (using https://data.london.gov.uk/dataset/land-area-and-population-density-ward-and-borough 2016)
library(tidyverse)
Density <- read_csv("E:/UCL/005-GI System&Science/GIS-Assessment 3/R assessment3/population_density_2016_final.csv",na = "n/a")
#normalization population_per_square_km
normalization<-function(x){
return((x-min(x))/(max(x)-min(x)))}
density.normalization <- normalization(Density$`population density`)
Density$Density.normalization <- density.normalization

#join the data to basemap P3P27
LondonWard.join.popdens <- LondonWardSF%>% left_join(Density,by=c("gss_code_ward"="New Code"))
qtm(LondonWard.join.popdens,fill="Density.normalization")

#### coexistence of old and new building
#read building list
building <- readOGR("E:/UCL/005-GI System&Science/GIS-Assessment 3/R assessment3/Listed Buildings/ListedBuildings_06Dec2018.shp",layer="ListedBuildings_06Dec2018")

library(maptools)
LondonBoundary<-unionSpatialPolygons(LondonWard,gss_code_borough)
plot(LondonBoundary)

POSAL <- read_csv("E:/UCL/005-GI System&Science/GIS-Assessment 3/R assessment3/access_public_open_space_ward.csv",na = "n/a")

Parsed with column specification:
cols(
  WD13CD = col_character(),
  `Ward name` = col_character(),
  `Borough name` = col_character(),
  `Open Space` = col_double(),
  `Local Parks` = col_double(),
  `District Parks` = col_double(),
  `Metropolitan Parks` = col_double(),
  `Regional Parks` = col_double()
)

POSAL$POSAL.normalization <- normalization(POSAL$`Open Space`)
LondonWard.join.POSAL <- LondonWardSF%>% left_join(POSAL,by=c("gss_code_ward"="WD13CD"))

Column `gss_code_ward`/`WD13CD` joining factor and character vector, coercing into character vector

#qtm(LondonWard.join.POSAL,fill="POSAL.normalization")
#combine two accessibility sub-index
Access.normalization<- merge(POSAL,PTAL, by.x="WD13CD",by.y="Ward Code",all=TRUE)

### Road intersections
node <- readOGR("E:/UCL/005-GI System&Science/GIS-Assessment 3/R assessment3/Road Networkoproad_essh_gb/new_road_node.shp",layer="new_road_node")
BNG = "+init=epsg:27700"
nodeBNG <- spTransform(node, BNG)
summary(node)
##qtm(nodeBNG)
nodeBNGSF <- st_as_sf(nodeBNG)
#select the node type"junction"
nodeBNGSF <- nodeBNGSF[which(nodeBNGSF$formOfNode=="junction"),]
nodeBNGSP <- as(nodeBNGSF, "Spatial")
head(node)
#select node point in borough
nodeBNGSP <- remove.duplicates(nodeBNGSP)
nodeBNGSP.final<- nodeBNGSP[LondonWard,]
# tmap view
tmap_mode("view")
tm_shape(LondonWard) +
  tm_polygons(col = NA, alpha = 0.5) +
tm_shape(nodeBNGSP.final) +
  tm_dots(col = "blue")

# count number of point in each ward 
library(GISTools)
poly.counts(nodeBNGSP.final, LondonWard) -> nodecount
##setNames(nodecount, LondonWard@data$nodecount)
London.node <- LondonWardSF
London.node$nodecount <- nodecount   #add number of point to each ward (London.node$nodecount)
#combine the ward area(sq km)
London.node$Square_Kilometres <- LondonWard.join.popdens$`Square Kilometres`
#calculate point/area--point per square kilometre
London.node$PointPerSK <-London.node$nodecount/London.node$Square_Kilometres

###economic vitality
#CLASS 1  house price
#HousePrice <- read_csv("E:/UCL/005-GI System&Science/GIS-Assessment 3/R assessment3/house_price_mean2016.csv",na = "n/a")
#HousePrice.normalization <- normalization(HousePrice$Value)
#HousePrice$HousePrice.normalization <- HousePrice.normalization
#CLASS 2 economic activity
EcoAct <- read_csv("E:/UCL/005-GI System&Science/GIS-Assessment 3/R assessment3/economic_activity_ward.csv",na = "n/a")
EcoAct.normalization <- normalization(EcoAct$`Economically active: Total`)
EcoAct$EcoAct.normalization <- EcoAct.normalization

recreationBNGSF <- union(recreationBNGSF2,recreationBNGSF,makeUniqueIDs = TRUE)

Error in union(recreationBNGSF2, recreationBNGSF, makeUniqueIDs = TRUE) : 
  unused argument (makeUniqueIDs = TRUE)

# Assign the land use area for each subset. Here we first use the most widespread residential as example 
library(raster)
r <- raster(ncols=400, nrows=400) #generate raster size(quantity)
extent(r) <- extent(residentialBNGSP)  #Important! Assign the extent of raster to cover the same extents of the polygon
residential.Raster <- rasterize(residentialBNGSP,r,background=NA) #convert polygon to raster
#tmap_mode("view")
#qtm(residential.Raster)+tm_shape(LondonWard) +tm_polygons(col = NA, alpha = 0.5)
#extract rasters from polygons#(extract() cannot use method=bilinear in raster-in-polygon)
residential.extract <- raster::extract(residential.Raster,LondonWard,df=TRUE, weights =FALSE, na.rm = TRUE)
#remove rows with NA
row.has.na <- apply(residential.extract, 1, function(x){any(is.na(x))})
residential.extract2 <-residential.extract[!row.has.na,]
#merge the pixel in same ward(aggregate the frequency of the same data in one)
library(plyr)
residential.extract.sum <- factor(residential.extract2$ID)  
residential.extract.sum <- table(residential.extract2$ID)
residential.extract.wardsum <- as.data.frame(residential.extract.sum)
#create new order number for LondonWardSF, for merge the residential.extract and LondonWardSF by order
LondonWard2SF <- LondonWardSF
LondonWard2SF$Var1 <- 1:630
LondonWard2SF[,'Var1']<-factor(LondonWard2SF[,'Var1'])
#merge two dataframe
residential.extract.final <- merge(residential.extract.wardsum, LondonWard2SF, by="Var1",all=TRUE)
#residential.extract.final [!duplicated(residential.extract.final ), ]
#replace NA values with 0
residential.extract.final[is.na(residential.extract.final)] <- 0
#rename the column
colnames(residential.extract.final)[2] <- "Freq.residential"
#raster cell size
#we have to convert raster to polygon because of raster(). This function only compute the pixel area in longitude/latitude coordiante system
residential.repolygon <- rasterToPolygons(residential.Raster, fun=NULL, n=4, na.rm=TRUE, digits=12, dissolve=TRUE)
#qtm(residential.repolygon)
residential.repolygonSF <- st_as_sf(residential.repolygon) 
residential.repolygonSF$area <- st_area(residential.repolygonSF)
#so the single pixel area is 28999.54 [m^2](0.029 km^2)

##based on the details above, we create functions for land use analysis
function.raster <- function(landsp){
  extent(r) <- extent(landsp)
  landsp.Raster <- rasterize(landsp,r,background=NA)
  return(landsp.Raster)
}

function.extract<- function(landsp.Raster){
  landsp.extract <- raster::extract(landsp.Raster,LondonWard,df=TRUE, weights =FALSE, na.rm = TRUE)
  row.has.na <- apply(landsp.extract, 1, function(x){any(is.na(x))})
  landsp.extract2 <-landsp.extract[!row.has.na,]
  landsp.extract.sum <- factor(landsp.extract2$ID)  
  landsp.extract.sum <- table(landsp.extract2$ID)
  landsp.extract.wardsum <- as.data.frame(landsp.extract.sum)
  landsp.extract.final <- merge(landsp.extract.wardsum, LondonWard2SF, by="Var1",all=TRUE)
  landsp.extract.final[is.na(landsp.extract.final)] <- 0
  return(landsp.extract.final)
}

qtm(residential.Raster)
library(tmap)
tmap_mode("view")
tm_shape(LondonWard) +
  tm_polygons(col = NA, alpha = 0.5) +
tm_shape(residentialBNGSP) +
  tm_polygons(col = "blue")